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Abstract 

A novel approach to unstructured quadrilateral mesh generation for planar domains is presented. 
Away from irregular vertices, the resulting meshes have the properties of nearly conformal grids. The 
technique is based on a theoretical relation between the present problem, and the inverse Poisson 
(IP) problem with point sources. An IP algorithm is described, which constructs a point-source 
distribution, whose sources correspond to the irregular vertices of the mesh. Both the background 
theory and the IP algorithm address the global nature of the mesh generation problem. The IP 
algorithm is incorporated in a complete mesh generation scheme, which also includes an algorithm 
for creating the final mesh. Example results are presented and discussed. 

1 Introduction 

Boundary alignment is a critical feature of meshes in many applications. In a boundary aligned mesh 
the boundary, or some other line, is traced by the sides of high-quality cells, see fig. [T] The definition of 
a "well-shaped" cell may be application dependent, but in many cases, cells similar in shape to squares 
(for quadrilateral cells) or equilateral triangles (for triangles) are preferred. Characteristics of the entire 
mesh are also important, such as smooth cells-size and cell-shape transitions. 

The problem of producing boundary aligned meshes with well-shaped cells has been the subject 
of extensive research [T]. Still, many popular algorithms are heuristic in nature, and a more general 
understanding of the subject is called for, especially when quadrilateral meshes are considered. A key 
difficulty is the problem's global character: the shape and position of every cell in the mesh is, at least 
in principle, related to that of any other cell. 

In a previous work [2] , we described a relation between the problem of two-dimensional unstructured 
mesh generation, on both planar and curved surfaces, and another well-known problem, namely the 
Inverse Poisson (IP) problem. The IP problem is concerned with reconstructing a source distribution 
p of the Poisson equation V^t/i = p, from information on the potential (j) at the boundaries. In that 
work, the mesh was assumed to be conformal away from the irregular vertices (vertices whose degree 
is different than four), like a grid mapped by a conformal mapping. Such grids have the property of 
having square cells in the limit of an increasingly finer grid. Under this assumption, the problem of mesh 



Figure 1 : Part of a boundary aligned mesh (black lines) . The gray line represents the domain's boundary. 
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generation was then shown to reduce to an IP problem. The irregular vertices of the mesh correspond 
to point sources (delta functions) of p, and (j) is interpreted as the logarithm of the local resolution. 
This theoretical framework turns the focus to the irregular vertices of the mesh: once their distribution 
is fixed, the continuum properties of the mesh - local resolution and directionality - are known. Note 
that it is also an explicitly global formulation, since the resolution at any given point is affected, via the 
function 0, by the locations of all irregular vertices in the mesh. 

In this paper the generation of planar quadrilateral meshes is discussed. Resting on the results of 
[2] , a new IP algorithm is presented, designed to construct source distributions of the appropriate type, 
which approximate the resolution and mesh directionality inputs at user-specified points, such as at 
the boundaries. The IP algorithm is then incorporated into a complete mesh generation scheme, which 
also includes a technique for generating the final mesh. An implementation is described, and shown in 
example cases to generate boundary aligned meshes, where well-placed point sources create smooth cell 
transitions and high quality cells. A similar procedure is probably applicable to triangular meshes, but 
is not discussed in the present work. 

Remeshing of curved surfaces has recently attracted considerable attention; for a review see [3] . Many 
of the algorithms receive an input mesh directionality throughout the surface, usually the principal 
curvature directions of the surface. This setting presents different challenges than those addressed here, 
since the mesh structure is determined, to a large extent, when mesh directionality is given everywhere in 
the domain. For example, the locations of critical irregular vertices are dictated by the mesh directionality 
in the vicinity of these point^. Another related subject is surface parameterization, concerned with 
creating mapping of surfaces to the plane. Conformal surface parameterizations are created in [5], but 
boundary alignment is not addressed. 

The paper is organized as follows: Section [2] shortly reviews the relevant theoretical background, 
with emphasis on the relation between the IP problem and unstructured mesh generation. Sections [3] 
and [4] describe the proposed IP algorithm, and the mesh generation scheme. Section [5] describes an 
implementation of the algorithm, and section [5] gives examples of meshes generated. Conclusions and 
possible directions for future research are discussed in section [7l 

2 Background theory, relation to the IP problem 

In this section an overview of the background theory will be given. Section 12.11 explains the rationale 
underlying the mathematical formulation. In section 12. 2i some key conformal geometry concepts are 
discussed, together with their relevance to the present problem. Section 12.31 summarizes the results 
developed in 2 , relating mesh generation with the IP problem. The exposition is limited to planar 
two-dimensional mesh generation; A detailed account, in the more general setting of curved surfaces, can 
be found in [2] . 

2.1 Motivation 

We start by considering mesh generation using conformal mappings, which can be viewed as a special case 
of the theory to be described. Mapping techniques, in general, construct a function from one domain, for 
which a mesh already exists, to a second domain, which is being meshed. The mapping function is then 
used to map the mesh into the second domain. A key idea is that continuum properties of the mapping 
function control the shapes of the cells of the new mesh, at least for small enough cells. For example, if 
the mapping is conformal, i.e. angle preserving, a cell with right inner angles (rectangle) will be mapped 
to a target cell with approximately right inner angles. 

In unstructured mesh generation the connectivity of the mesh is not known in advance, and a more 
general framework is called for. In what follows, the interplay between the two domains which serves as 
a paradigm of mapping techniques, is replaced with an interplay between two definitions of distances on 
the input domain. Instead of imagining the mesh as being the image of a mesh on a different domain, 

^In many works, both mesh directionality and local resolution are specified on the entire surface, and since conformality 
requires a specific relation between the two properties, the resulting meshes are not conformal (nor are they claimed to 
be). In one exception [4], a preprocessing step attempts to adjust the local resolution, in order to create a more closely 
conformal mesh. 
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Figure 2: a. A conformal grid. Solid lines are geodesies. The distance between two parallel lines is ho, 
measured with the new distance definition. Geodesies bend towards smaller cells. Crosses represent the 
cross-field at selected points, b. The curvature of a geodesic is related to the gradient of 0. 



we redefine the distances on the domain to be meshed, and fix the cell edge length. Thus, the local mesh 
resolution will be proportional to the new local distance definition: a large new distance between two 
given points will mean more cells will be placed in between, hence a higher local resolution. Distances 
are redefined using the concept of a metric, known from Riemannian geometry. 

Since we focus on cells which are squares in the limit of an increasingly finer mesh, we need to define 
only two local properties: the resolution (inverse of cell-size), and the direction of the cell-edges. The 
new resolution, as noted before, is controlled by a new distance definition. We use a new distance which 
is localy a scaling of the old distances: the new distance between a point and other nearby points is equal 
to the old distance, multiplied by a scalar factor which is independent of the direction. Thus, a small 
square measured in one distance definition is also (approximately) a square according to the new distance 
definition. In Riemannian geometry terminology, the old and new metrics are said to be conformally 
related. The mesh directionality is related to the resolution, as is described in the following section. 

2.2 Definitions, conformal geometry relations 

It is convenient to work with a function (/) (r), defined as the logarithm of the local scaling factor. That 
is, a small square of side length h measured with the original, Cartesian distance definition, is a square 
with side length h{r) = /i (r) e"^'""' as measured with the new distance definition, see fig. Oa.. If we 
imagine that the domain is covered with many small squares, all with the same side length ho in the new 
distance definition, the size of these cells in the original distance definition will be h (r) = hoe~''^^^\ The 
(local) resolution is the inverse of the local-cell-size, hence 

resolution (X e'^'-^^\ (1) 

where the proportionality constant is a single number for the whole mesh. 

The second continuum property is the local directionality of the cell edges. The edges should ideally 
meet at right angles everywhere, except at irregular vertices, where the number of edges incident on the 
vertex is different than four. It is therefore natural to assign to every point a set of four directions, 
mutually parallel or perpendicular. This concept was expressed by many authors, and given various 
names, such as mesh directionality [6], 4-symmetry direction field [7], and frame field, although the 
last may refer to a structure which also holds cell-size information [5]. Graphically, this object can be 
represented by a cross at every point, see fig. [Ha., and here will be called a cross-field [2J. On the plane, 
the cross direction can be measured by the angle 9 from the x-axis to one of the directions of the cross. 
This angle is fixed up to an addition of tt/2 radians, i.e. 9i, 62 represent the same cross iff 9i = 92+mT/2, 
for some integer n. 

The function ip and the cross-field are not unrelated; due to conformality, lines that trace the edge 
directions bend towards the side with smaller cells, see fig. [21a.. In the continuum theory, these lines 
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are known as geodesies of the manifold. Geodesies are a generalization of the concept of straight lines 
to non-Euclidian geometries. In the original, Cartesian coordinates, the geodesic obeys the following 
differential equation: 

where k is the curvature of the geodesic in the Cartesian coordinate system, and d4>/dn is the derivative 
of the function cj) in the direction normal to the tangent, see fig. Ob.. Eq. ([2]) allows one to calculate 
the change in the direction of the cross-field between two points connected by a geodesic. There is also 
a direct way of calculating the change in the direction of a cross-field between any two points, along any 
curve connecting the two points, known as parallel-transport. Let a be some curve from point a to point 
b, and let 9a, Oi, be the angles of the crosses at points a,b respectively (as noted above, the angles are 
defined up to an addition of a 7r/2 radians). Then 

^ds = 9t~9a. (3) 
on 

where the integration denotes a line integral along the curve a, according to the length parameter s on 
a, as measured in the original coordinate system. The differential formulation of this equation is 

^ = ^ 

dn dt ' ^ ' 

where (t, nj arc a pair of perpendicular vectors that form a right hand system. 

Where the function is defined, that is, at any point in the domain that is not a singularity, cj) can 
be shown to harmonic, that is to obey the Laplace equation: 

where, again, the derivatives are taken with the original coordinate system. Eq. (|2l5p are well-known 
results of conformal geometry [9], [10]. Eq. ([3]) is derived in ^ 



2.3 Relation to the IP problem 

Equations (|2l5p fully describe the relations between the cross-field and the function <j), at any regular 
point of the domain, that is, any point that is not a singularity of the function cj). The singularities of the 
function (j) are a key ingredient of the theory, since they correspond to the irregular vertices of the mesh, 
and unstructured meshes are those which contain irregular vertices. A detailed analysis of the possible 
of singularities of the harmonic function (jj, and their effect on the resulting mesh, was carried out in [2], 
and shows that the only type of singularity that corresponds to a mesh with a finite number of cells is 
of the type (f) (r) oc In |r — ro|, where tq is the location of the singularity. In the IP literature, such a 
sigularity is known as a point source. Furthermore, the prefactor of the logarithm is directly related to 
the degree of the irregular vertex in the final mesh. More specifically, suppose there are ric singularities 
in a domain D, at points r,„, m — L.Uc, then the function (j) can be written as 
Condition 1: 

71(2 

0(r) = ^ In |r - r„i\ ^ (t^L + -^/cmln|f - frn\ , (6) 

m— 1 m=l 

where (f>L is a harmonic function, km S and km > —4. The numbers Qm are known as charges, so 
Condition 1 states that the charges are integer multiples of tt/2. We will refer to the numbers km simply 
as the "fc -values" of the singularities. The degree (number of incident edges) of the irregular vertex 
corresponding to the source at is equal to 4 -I- km- So, for example, irregular vertices of degrees 3 and 
5 will correspond to a singularities with km = — 1 and km = +1, respectively. 

In a geometrical context, a logarithmic singularity of (j) represents a cone-point - the tip of a cone - in 
the new distance definition. The charge corresponds to the angle deficit of the tip of the cone. Related 
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a. Condition 3: b. Condition 4: 





Figure 3: Conditions 3 and 4. 



subjects where manifolds with cone-points are considered include the study of disclinations in elastic 
media [TT], and surface parameterizations [T^ . 

An unstructured mesh is required to be aligned with the boundary. Namely, that one of the cross 
directions be parallel to the tangent to the boundary (where a tangent exists) : 

Boundary Alignment definition: 



tan (e (p) + n|) - ^ (p) , 



per. 



(7) 



where (y (s) , a; (s)) is a curve tracing the boundary, p = {x,y) and n e Z. Since 9 is fixed up to an 
addition of a multiple of 7r/2 radians, this definition depends only on the cross-field itself, and not on 
the particular choice of 6. 

Requiring boundary alignment, as defined in eq. ([7]) can be shown to be equivalent to three more 
conditions on the function [2, . The following. Condition 2, applies to a point on a smooth section of 
the boundary, r. It roughly states that a smooth section of the boundary is a geodesic. This provides 
information on the derivative of normal to the boundary: 

Condition 2: 

^{p)^K{p), per. (8) 

K is the boundary curvature at p. The next condition is concerned with junction points of the boundary, 
where the tangent to the boundary is discontinuous. It assures that the cross-field will be aligned with 
the boundary on both sides of the junction point. Let a be a curve from point a to point b on both sides 
of a point c of the boundary, see fig. Ob.. Then 
Condition 3: 

— ds = n- + 6'i„, (9) 
_ _ on 2 

where 6in is the angle of the boundary, and n G Z. This means that for some inner angles will contain 
a singularity at c; at a distance r from c, the singularity will be of type d(j)/dr ^ 1/r. 

The final condition is concerned with the relation between two different boundary components. It is 
just a restatement of eq. ^ for two points on two different boundary segments. 

Condition 4: 

Let a, 6 be two points on two boundary curves r(j,r^. Let o. be a curve from a to h. Xhen 



— ds 
on 



(10) 



Conditions 1-4 form a complete set of conditions on 0, for a cross-field that is boundary aligned to exist. 
That is, conditions 1-3, and condition 4 on a number of selected curves, one from a selected point to 
each boundary component, are sufficient for a boundary aligned cross-field to exist, see [2]. 
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3 Algorithm 



3.1 Mesh Generation Steps 

The mesh generation scheme consists of the following steps: 

1. Setup and solution of the IP problem. Finds the sources' locations and charges. Described in 
sections 13.2113.61 

2. Solution of the direct Poisson problem, to obtain the functions (p, 9 throughout the domain, see 
section 13.71 

3. Generation of the final mesh, see section IH 



Suppose that we are given an open domain D, with boundary F, and are given the required resolution 
function on F. Using eq. ([T]), this is readily translated to the required value of on the boundary: 



Condition 1 (eq. ([6])) is a solution of the Poisson equation V^i/) — p, with point sources, whose 
locations and charges are yet unknown. The problem is to find a distribution of sources (location and 
charge) adhering to the boundary alignment definition, eq. ([7]), or alternatively to Conditions 2-4 (eq. 
([8|)- ([T0|) ). as well as to eq. (fTTj). Such a problem is known as an Inverse Poisson (IP) problem. The 
IP problem may be compared with the direct Poisson problem, where the source distribution is given, 
as well as some boundary information (e.g., Dirichlet or Newmann boundary conditions), and the value 
of the function (j) is to be found. In the IP problem, the source distribution is unknown, and a source 
distribution adhering to the known boundary information is to be found. 

IP problems have important applications in various areas of science and engineering |13j.|14j.[15j.[16j 
■ |17|. By its nature, the IP problem is ill-posed, and the solution might not be unique, and may be 
sensitive to small changes of the input, such as small changes in boundary conditions. In delicate problems 
of this type, any prior information on the source distribution may greatly affect the applicability of a 
specific solution procedure. In the present problem, we seek a point source distribution. A number of 
algorithms for solving an IP problem with point sources appear in the literature. In |18j . an inverse 
problem where all sources have the same charge is solved. In our implementation, however, at least two 
different charge values must be incorporated {Q = ±7r/2). In [IS] , [10] , [H] an inverse problem where both 
the locations of sources and their charges are unknown, and are reconstructed. This gives more freedom 
in the reconstruction than we can allow, since for the present purposes the charges must be multiples 
of 7r/2. Another important aspect of the present problem is that the domain of the IP problem may be 
of any shape and topology (i.e., may contain holes), whereas the above works only deal with a simply 
connected region (a circle, usually). 

It is important to note that in the present application, in contrast to other standard applications, the 
input to the IP problem is not generated by some existing source distribution (perhaps with some added 
noise), but by the domain's shape and input resolution. The existence of a source distribution which 
reconstructs the input data, at least approximately, is therefore not obvious. This is an interesting and 
important subject, but is beyond the scope of the present work. 

3.3 Complex Formulation 

As is well known, the real and imaginary parts of a complex analytic function are harmonic functions 
|22j . This correspondence has been utilized in IP algorithms [20] , [19] , [21] , and will be used here as well. 
We define the complex-valued function 



3.2 



Input 



In (resolution)\Y, . 



(11) 




(12) 



m— 1 
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where z,„ = x„i + iym, with {xm,ym) the components of fm, that were defined in Condition 1. /i is a 
function on D, such that 

Re (h) — (j)L- 

Then, recaUing that Re (In (z)) — In \z\, it foUows that 

Re(F(z)) = 0(f), (13) 

with (j) (r) defined in Condition 1, and z — x + iy, with r = {x, y). 

The functions h{z) and In (z — z,„) for some Zm are analytic in a neighborhood of any point that 
is not a singularity of cj). However, as functions over the entire domain D, they may be multi-valued. 
Multi-valued functions accept many-values at a point, depending on the path taken to that point. It is 
well known that the complex function In (z) is multi-valued. Defining In (z) as ln(z) = J^dt/t, where 
7 is a path from 1 to z, the imaginary part of ln(z) is only fixed up to an addition of multiples of 2tt, 
depending on the path taken. If the domain D is not simply-connected (i.e., contains holes), then the 
function h (z) may also be multi-valued, since by following different paths around the holes, different 
values of h (z) may be obtained. 

The function F (z), being a sum of multi-valued functions, may also be multi-valued. According to 
eq. (|13p . the real part of F is single- valued, as it is equal to (f> at that point. We thus turn to examine 
the imaginary part of F. The real and imaginary parts of a complex function are related by the Cauchy- 
Reimann (CR) equations [52]. Recall that for a complex differentiable function / (z) = u{z) + iv{z), 
where u(z),v{z) are real-valued, the Cauchy-Reimann equations read d^u — dyV, and dyU = —d^v. 
This allows one to recover the imaginary part of a analytic function if its real part is given. Writing 
eq. ([ll in the two right hand coordinate systems {x,y) , {y, —x), gives the relations d(j)/dy — dO/dx and 
d(f>/dx — —dO/dy. These are precisely the CR equations for the complex function </) — i (6* + C), where 
C is any real constant, hence 

F{z) = 4>{z)-i f ^ds = 4>{z)~i{e (z) + C) . (14) 
J~f on 

The constant C is arbitrary, and will be taken to be zero. The multi-valued nature of F is explicit in 
the integral formulation in eq. (HH). The right-hand-side (RHS) of eq. shows that the since the 

imaginary part of F is —0 (z), and 9 must be fixed up to additions of 7r/2, then the integral 

must also be fixed up to additions of n/2. Indeed, this was shown to hold if Conditions 1-3 (eq. ©-([H])) 
hold, see [2]. 

To summarize, the problem is restated as finding a multi-valued complex function _F, given by eq. 
(fT2|) . and whose value at points on the boundary of D is: 

F{z)^(j){z)-ie{z). (16) 

In the following sections, an algorithm for solving the problem as it was here restated is described. 
The IP algorithm first removes the contribution of h (z) from the boundary conditions, paying special 
attention to the junction points, and then constructs the source distribution. 

3.4 Handling junction points 

The function and hence the function F ^ may be singular at junction points. Special procedures for 
addressing this behavior are taken as part of two different steps of the algorithm: 

a. Choosing the input resolution near junction points. 

b. Special treatment of junction points when h is removed from the boudary conditions. 
The two subjects are discussed in the following subsections. 13.4. II and [51^1^ 
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3.4.1 Choosing the input resolution near junction points 



Except in the special cases when 9 in is a multiple of 7r/2, the behavior of at a small neighborhood 
in D of a junction point is singular, with the singularity at the junction point, see Condition 3 (eq. 
([5])). If the input resolution in a small neighbourhood of the junction point does not match this singular 
behavior, solution of the IP problem will feature sources at any distance from the junction, no matter 
how small, resulting in a distribution with an infinite number of sources. In practice, this means that 
an IP algorithm will cluster many sources near the junction point, in a futile attempt to reconstruct 
the boundary conditions there. To avoid this problem, we adjust the input at a neighborhood of the 
junction point. This adjustment should vanish rapidly at a distance larger than 1-2 cell-sizes, so as to 
have little effect on the original required resolution. 

Condition 3 states that the gradient of near a junction point must diverge as d(j)/dn ~ 1/r, where r 
is the distance from the junction point. Such a fiux is formed by a logarithmic singularity at the junction 
point (see also section 7.1). Consider a singular source term of the form 4>{r) = ^Inlrj, where r 
is the distance from the junction point. Then the flux through a circular arc a^j at distance r from the 
junction point is (see also fig. [3la.): 

on ZTT 

where 9in is the junction inner angle. Then by Condition 3, = fc-| -t- 9in, hence 

Qj = 2n (nj^ - l) . (17) 



nj is integer, which must be positive, otherwise a mesh with an infinite number of cells will result ([5], 
section 7.1). In fact, nj has a simple interpretation: it is the number of cells incident upon the junction 
point in the resulting mes A reasonable choice for nj is therefore 



nj = round j > (1^) 

for which the inner angles of the cells incident on the junction point are closest to 7r/2. 

In order to restrict the effect of this correction of to a small region in D, a source term with the 
opposite charge —Qj can be placed outside D, at a distance of about 1-2 edge lengths. See for example 
the second and third examples in section [6l 



3.4.2 Treatment of junctions when removing the harmonic part 

In section [3.51 a method will be described for calculating and removing the contribution of h (z) to the 
boundary conditions. At junction points, however, the function h (z) may be singular: as D is an open 
set, boundary singularities are part of h{z), so care should be taken when preforming the calculations 
described in section [3.51 especially in a numerical implementation of the technique. 

To avoid these problems, we subtract the junction singularities from the boundary value of F [z) before 
proceeding with removing h [z). Let z^ be the locations of the junction points, and Qj. the charges as 
given by eq. (1171 ). (jlSp . The complex function corresponding to In \z\ is In (z), so the contribution of the 
junctions is 

-EI^^I- (19) 

If a singularity is located on an inner boundary component, the logarithms must contain a branch cut 
somewhere in D. In order to avoid adding a cut, we add another term to every hole, equal to minus the 

^During the creation of the final mesh connectivity (section|4]l nearby cone-points may be joined, and cone-points may 
be shifted to the boundary, depending on the final cell-size. In such a case, the number of cells incident on a junction point, 
as well as on other boundary points, may change. 
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Figure 4: Removing the harmonic part using the Cauchy integral, a. The case of a simply-connected 
domain, b. A domain with more than one boundary component. Thick line represents the boundary, 
thin line the cut introduced. 



total charge in each hole. Denote the inner boundaries by .., r„^^^, and let z!^"'-^, . 
points inside the respective holes. The following terms are therefore added to F (z): 



^hole 



be arbitrary 



Fiz)~.F{z) 



— — In (z — z 



In (. 



hole\ 

Z Z^ j 



(20) 



m=2 JiGT. 



In this way, the additional terms form an analytic function in D, with no branch cut, and are a part of 
h{z). Since the contribution of h{z) is removed from F (z) as described in the next section, the exact 
form of this additional term, i.e. the choice of points z'^''^, does not affect the results. 

3.5 Removing the harmonic part 

In this section a technique for calculating h (z) is presented. After h (z) is calculated, its contribution to 
F (z) at the boundary can be subtracted. To simplify notation, let Fp (z) denote the sum of logarithms 



In the context of IP algorithms, the idea that the contribution of the harmonic part to the solution 
can be removed from the boundary information was suggested in [21j. There, an IP problem in the unit 
disk B is discussed. The value of the analytic part of F (z) on the boundary of B can then be calculated 
by taking the Fourier transform of F (z), and leaving only the positive frequency components, as can be 
shown by considering the Lourent series of F (z) in the unit disk. By its construction, this technique 
applies to functions defined on the unit disk. It can be extended to other domains, if a conformal 
mapping of the domain to the unit disk (which exists according to the Riemann mapping theorem) is 
calculated. For our present purposes, since the domain D in every problem input is different, using 
this technique would require constructing a mapping to the unit disk for each domain being meshed 
separately. Furthermore, the domains in our problem may contain holes, which further complicates the 
matter. We therefore use an alternative approach, which is now described. It is based on the Cauchy 
integral theorem, and is similar to the "sum factorization" step in the Weiner-Hopf technique [53] , applied 
to a bounded domain (O, defined below) in place of an infinite strip. 

As in [21] , we assume there is a positive distance d between the boundary of D and the source closest 
to it, as is the case, for example, when the number of sources is finite. Let be the set of points at a 
distance smaller than d to the boundary, see fig. (Ha.. As will become apparent later, the value of d does 
not enter the calculation and is irrelevant, as long as there exists some d > as required. 

We first consider a simply-connected domain D, with boundary F. Define the fiux through any curve 



in F(z): 




(21) 



m— 1 



a: 




(22) 



9 



The flux $r through the boundary can be calculated: the value of d(f>/dn at smooth boundary points is 
known from Condition 2 (eq. ([8])), and the flux at junction points is calculated from Condition 3, as is 
explained in section [3.4. II above. This flux can be shown to be a multiple of tt/2. 

Remark 1 Using Conditions 2,3 and the fact that the rotation of a tangent in a simple curve is 2tt , the 
total flux through a boundary component r„ can he shown to he equal to s^(— 4 — X]j er ^''^Ji ~ ^))' 
where Ji are the junction points, uj. for each junction is given hy eq. I118\} . and s — +1, —1 for inner 
and outer boundary components, respectively. See also f^. 



To proceed with the decomposition calculation, we assume that $r = 0, as calculated from eq 
above. If this is not the case, following [21] , we subtract a source term to the boundary conditions, 
centered at some point Zin inside the domain: 

F^F-^ln(z-z„) (23) 

such that $r is zero . As will be shown below, the choice of Zi„ does not affect the results. Since $r 
was equal to a multiple of 7r/2 before, kin G ^- When <i>r = 0, the function F (z) is single-valued and 
analytic in SI, as follows from the integral representation of F (z) in eq. p^ . along with <i>r = 0. We 
now use the Cauchy's theorem, stating that for every point a £ SI: 

F{a) = ^ I ^dz. (24) 
' 2ni Jg^z-a ^ ' 

The boundary dU, of SI has two boundary connectivity components: the outer boundary (which is also 
the boundary of D), and the inner boundary. Denote them by d^lo and dVli respectively, see fig. |31a.. 
Eq. ((24l) now reads 

Fia)^^! ^dz + ^[ l^dz. (25) 
27r« Jdn^ z-a 2m Jg^^ z-a 

The function Re (F) — (p is harmonic on SI. According to a decomposition theorem (see [H], Chapter 9) 
its decomposition into Re (Fp) and Re (h) — (f>L, is also unique. It follows that the decomposition of the 
multi- valued complex function F into Fp and h is unique, up to the arbitrary constant C in eq. (|14p . 
chosen before to be zero. The two integrals expressions on the RHS of eq. ((25|) correspond exactly to 
the two components of the (unique) decomposition of F (z) in Q as described in |,24j , hence 



h{z) = ^. I ^dz; Fp {z) = ^ I ^dz, 
Zm Jy z ~ a Zm J qq^. z — a 



(26) 



where we have used the fact that F = 9Slo. We would like to find Fp on F, but dQ.i is not known, nor is 
F {z) on (9S7i, so Fp (z) cannot be computed directly from the second equation. The first equation can 
however be used, since F (z) is given on F, so h (z) can be calculated on F (more precisely, since h is 
defined in Z3, it is the limit of /i as F is approached). Then Fp (z) on F is given by Fp (z) — F (z) — h (z), 
according to eq. (pij) . 

Finally, we add back the source term subtracted before 

Fp ^ Fp + ^ In (z - z„0 . (27) 

We now turn to the case when D is not simply-connected. We add source terms: one inside the 
domain, as discussed above (eq. ([^5]) ). and one inside every hole (i.e. outside D), such that the fiux $r„ 
through every boundary connectivity element is zero, see Fig. [Hb.. 

Fp^Fp-J2^^^{-- 4:'^) (28) 

m=2 

The z'^'-'^ can be the same points used in section [3.4.21 or other points. The results do not depend on the 
additional charges' locations, as will be explained below. We now introduce cuts so that the boundary 
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contains a single connectivity element (the cuts play a part in the derivation, but drop out of the final 
calculation), see fig. [Ub.. The cuts introduced also serve as the branch cuts of the logarithms. Denote 
this new boundary F. Using F, we proceed as when D is simply connected: define Q as before, see fig. 
lUb., and subtract the source term inside D, as in eq. ([23|l . According to eq. ([26]) : 



,/N 1 f F(z) , 1 f F(z) , 1 f F(z) , 

'^i^) = 7;- L-^dz^— -^dz + — -^dz. (29) 



2TTi Jf z ~~ a 2711 Jy z — a 27ri 

The second integral on the RHS denotes the integral over the cuts introduced to form F. Note that each 
cut-path is traversed twice, back and forth. Since the flux through each and every hole boundary 
is zero. The value of F {z) ^ (j) — i J dn(f>ds, when integrated along F, is continuous across the branch 
cuts, and the integrations over each cut traversed in both directions cancel each other, and drop from 
the total integration. Therefore 

h{z) = ^ f ^^dz (30) 
^ ' 2m Jy z - a ^ ' 

as in eq. ([26|) . As before, Fp (z) — F (z) — h (z). Finally, the source term inside D is added back, as in 
eq. ([27|) . The boundary value of Fp obtained after this source term is added back may be multi- valued, 
with a branch cut discontinuity in the imaginary part that is a multiple of tt/2. This is a valid input to 
the next step, as explained in section [3T6l below (following eq. ((3T|) '). 

The choice of the added singularities' locations outside D (inside the holes) does not affect the results: 
changing the location of a singularity with some charge from Za to Zb is equivalent to adding a singularity 
with opposite charge at Za, and with the same charge at Zb- Since these two singularities lay outside 
fi, and are of opposite charge, this amounts to adding an analytic function to F, which is removed by 
the Cauchy integral technique described above. The result is also unaffected by the location Zi„ of the 
source term subtracted inside D, since its location only affects the Fp and it is later added back. The 
choice of cuts in eq. [29l of course, does not affect the Cauchy integral calculation, since the cuts do not 
enter the final calculation, eq. ((30|) . 

3.6 Calculating the sources' locations 

At this point, we are given the value of Fp (z) on the boundary F of the domain D. Using eq. (fT2|l . we 
define 



^ (z) = exp [AFp (z)] = exp 



km In {z 



,m—l 



(31) 



Note that while Fp may be have a multi-valued imaginary part with a branch cut discontinuity: i6 —^ 
iO + m7r/2, for some integer n, ^ (z) is single valued, since exp [4 (m7r/2)] = 1. 

Since the fc-values, fcm, are integers, we can assume without loss of generality that km = ±1. Other 
charges may be formed by placing several sources at the same location. We can therefore write eq. \12\ 
as: 

AFp (z) = In (z - z+) - ^ In (z - z„) , (32) 

m—l m— 1 

where n^-^n^ are the number of sources with km = +1,^1 respectively. 

Remark 2 According to the divergence theorem, the difference — n_ is equal to $r/ (""/S) (this can 
also be shown directly by calculating the integral of Fp along T). Since $r is known (see also remark 
QP, — n_ is known. Therefore, to fix n^,n^ only one number remains to be chosen, e.g. n_|_ -|- n_. 
The choice of the total number of sources affects the accuracy of the reconstruction of the input data, see 
sections \3.8\ and section\^ 
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Substituting eq. ([32]), eq. ((3T|) now reads 



n (^-^™) 

ii^)-"^ • (33) 

n - 

m— 1 

In eq. (j33p . the numbers are the roots of polynomials, and we proceed by presenting the 

polynomials in a different form. Later, once the polynomials are found, the are recovered by 

finding the roots of the polynomials. We rewrite eq. (|33p as 

71+ —1 

= ' (34) 

ni—0 

where pi, qi are unknown coefficients. Note that the prefactors of z"+, z"" are indeed equal to 1, as can 
be seen by expanding the polynomials in eq. (|33p . To find the unknown coefficients numerically, we 
discretize eq. (p4l) . and give ^ (z) at N points Zj, j = l..iV, and write = ^ (zj): 

— 1 

The RHS of eq. (|35p is a rational Junction interpolation of the data (see e.g. 25 ). In the im- 
plementation described below, the unknowns Pm,<lm are evaluated as follows. Rearranging, eq. (|35p 
reads 

— 1 n_— 1 

?j ^ 2^P^ ) ~ 1^1^ [^3^] ) ■ (36) 

m— m— 

These are TV linear equations for the + n_ unknowns: Pm, Qm- Since TV will typically be larger than 
7i_i- + n_, the solution will only be approximate, e.g., a solution in the least- mean-square (LMS) sense. 

Once the unknown variables Pm,<lm are found, the source terms' locations z~,z^ are calculated by 
finding the roots of the two polynomials appearing in the RHS of eq. psp . 

Remark 3 A solution in the LMS sense, as described in eq. i36]) . is not recommended if the resulting 
errors in the 's are large, when compared with . This is due to two reasons: First, the input data 
appears in both sides of eq. i36\) . so the errors in the LMS solution might not reflect the errors in ^j, and 
moreover, the LMS solution may be sensitive to the choice of the coordinate system origin. Secondly, an 
error in scales as a difference in the resolution (since \ = exp((j})), rather than the more "natural" 
error definition given by the ratio of input and obtained resolution. This means that at low resolutions 
(large cell-sizes), errors in ^ might represent large relative resolution deviations. 

3.7 Restoring the harmonic part 

Once the sources' locations have been determined, the function (f> can be calculated everywhere in the 
domain. The solution will usually be approximate, i.e., the resulting sum of sources will only approximate 
the required Fp = F — h, and some deviation from both boundary alignment and resolution requirements 
may be found. The trade-off between the two requirements can be partly controlled by choosing an 
appropriate harmonic 0^ in eq. ©. In order to satisfy the cell-size requirements exactly, one can solve 
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the (direct) Poisson equation, given by eq. ([6]), but this time with a known charge distribution, and with 
Dirichlet boundary conditions, i.e. by specifying cj) on the boundary. 

A different approach is to try and satisfy the boundary ahgnment conditions, by solving eq. ^ with 
Newmann boundary conditions. This gives a new on the boundary. The difference between the input 
(required) </), and the (/> obtained on the boundary can be used to evaluate the quality of the source 
distribution calculated by the IP algorithm, and the total number n+ + n_ of sources can be changed 
accordingly. Note that when D is not simply connected, solving 4>l with Neumann boundary conditions 
(Condition 2) does not mean that exact boundary alignment is obtained, because Condition 4 might not 
be exactly fulfiled. 

3.8 Summary of the IP algorithm 

The steps followed in calculating the source distribution (locations and charges) can be summarized as 
follows: 

1. Adjusting the input resolution near the junction points fsection l3.4.1[) . 

2. Subtracting the junction source terms (sections 13.4.213.5)1 . 

{ju7ictions^ m—2 Ji^Tm 

-^ln(z-z,„)-^^ln(z-zri- 

Tn=2 

3. Calculating the Cauchy integral to find h (z) on the boundary. Finding Fp = F~h on the boundary, 
(section 13. 5( ). 

4. Adding back the source inside D (section [3. 5 1) : 

Fp^ Fp + ^ln{z- Zin) ■ 

5. Calculating the locations of sources fsection l?^ , using some initial total number of sources n+ +n_ . 

6. Solving the (direct) Poisson problem to find the new F approximating the original F. According 
to the quality of reconstruction, step 5 can be repeated with a different total number of charges. 

4 Creating the final mesh 

Once the function (j) is set throughout the domain, the final mesh can be constructed. There are various 
possible approaches to this problem. We present a simple method that was used to create the examples in 
section |6l In this method, the domain is cut along geodesies that follow the cross-field directions. Every 
source is at the end of a cut, and every inner boundary component is connected to the outer boundary 
by a cut, see fig. [Ua.. The direction of the geodesies emanating from the sources and directed along 
the cross-field are calculated as explained in Appendix A. The cut domain, which we denote by D' , does 
not contain any sources in its interior, and its boundary consists of a single connected component. A 
conformal mapping g of D' into the plane can be constructed, with a conformal factor (local scaling) of 
exp(0). 

We consider D' and g {D') as lying in the complex plane. First, note that since D' contains no 
sources in its interior and is simply-connected, (z) is uniquely defined in D'. Denote by 9d' (z) this 
single- valued 9 (z). To calculate the mapped boundary g (£>'), we calculate the function g (z) along the 
(single) boundary F' of D': 

g {z)\-p, — J exp {ip {z) ~ iOo' {z)) dz. 
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a. b. 




Figure 5: Creating the final mesh. a. The cut-tree. Points a, c are source locations, b. The mapped 
cut-tree. c. Modifying the boundary. Integer grid lines are drawn in gray. Sections a6, be, cd, de and ef 
■ATc shifted to the nearest grid line (dashed line) . Section fg is shifted so that de and ef will be of equal 
length, d. Section gh is shifted so that the lengths bc = cd+ fg. The process is continued until point a 
is reached again. 



14 



If this mapping is used to create a mesh, by placing a Cartesian grid on g {D'), and mapping it back 
to D' by taking the inverse of g, the mesh obtained wiU be invahd: the grid wiU not be continuous across 
the cuts, will contain invalid cells at the sources, and cut cells at the boundaries. To correct this, and 
obtain a valid final mesh, we seek a function / (z) on T' approximating g (V). A related problem, for 
surface meshes without boundary alignment, was addressed in [5l. 

First, note that because the cuts are added along geodesies aligned with the cross field, g (F') is a 
piecewise straight line, whose straight segments are either horizontal or vertical. In order to obtain a 
valid final mesh, the new function / (z) should have the following properties: 

1. / (z) must map the closed path F' to a closed path. 

2. / (F') must be a piecewise straight line, composed of horizontal and vertical lines only. To each 
segment in g (F') corresponds a segment in / (F'), at the same order along the paths. 

3. Each segment in / (F') must be directed in the same direction as the corresponding segment in 
g (F'), or be of zero length. 

4. If 7' is a geodesic cut of F', it is followed in both directions: 7' and 7'^. Then / (7') and / (7'^) 
must have the same length. This ensures the same number of cells on both sides of a single cut. 

We describe a simple algorithm for creating / (F'), meeting these requirements. To create / (F'), the 
path (? (F') is followed, and the segments of g(z)|p, are modified one after the other, in the order they 
appear in g (F'), see fig. [SJ Each segment is shifted to a parallel line of the integer grid. This is usually 
the nearest segment, except in the case the length of the segment directly before this segment dictates a 
different shift, according to properties (3,4) above, see fig. Elc.. This algorithm requires that at least one 
junction exist in F, which is mapped to a right angle in g (F'). The algorithm starts from the segment 
following this junction, so that when the last segment is reached and shifted, a change in the length of 
the first segment will be allowed. 

The resulting modified path is similar to g (F') , only with different segments lengths. / (F') is 
then defined as the composition of g (F') followed by a linear mapping of each segment of g (F') to 
the corresponding segment in the modified path. Once / (F') is known, / (£)') can be defined. In the 
implementation described below this is done by placing a triangular mesh in D' , and solving the Laplace 
equation for both Re (/ (z)) and Im (/ (z)) on D' . f (z) obtained in this way is not in general analytic, 
since the CR conditions might not hold between Re(/(z)) and Im(/(z)), so the final mesh will not 
exactly conformal, but only approximately, as f (z) approximates the analytical function g{z). The 
mesh edges are then extracted by tracing the integer- valued lines of Re (/ (z)) and Im (/ (z)). 

Once the connectivity of the mesh has been established, a final smoothing procedure which involves 
all the interior vertices, including those on the cut path, is preformed. In the definition of / (z), both 
Re(/(z)) and Im(/(z)) are harmonic functions. A standard finite difference approximation to a har- 
monic function on a uniform grid, assigns to each vertex the average of its four neighboring vertices. 
This is applied to all vertices inside the mesh, including those on the cuts of F', which can be interpreted 
as a generalization of the defintion of / (z) (to be more precise, in this smoothing procedure it is (z) 
which is assumed to be harmonic). This procedure is the well-known Laplacian smoothing procedure, in 
which the location of every vertex is equal to the average location of its neighboring vertices, and in this 
application it is theoretically justified if / (z) is approximately g (z), i.e. if / (F') is close to g (F'). 

5 Implementation Details 

The algorithm described in the previous section was implemented in Matlab [27 , using standard Matlab 
functions, such as a least-mean square solver for linear systems, an ODE solver, etc.. 

To calculate the Cauchy integral (section 13. 5p , the boundary was approximated by a piecewise linear 
path between sampled points, with a constant value of 4> on each segment. This allows the integral over 
each segment of the path to be calculated analytically, and the path integral there is equal to the sum 
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of the integral over the segments. This is a low order approximation, but was sufficient for the present 
purposes. 

The linear system, eq. ([36|). was solved in the least- mean-square sense, using Matlab's mldivide 
function. Note that the coefficients in this system contain powers of the Zj's, the boundary data locations. 
This leads to ill-conditioned linear systems, which are more sensitive to the numerical round-off errors, 
as the number of sources n_|.,n_ increases, or the range of |z| increases. Tests indicate that when the 
solver issues a warning, indicating that the matrix is rank deficient to the working precision, the results 
of the IP algorithm become unreliable (this was not the case in the examples in section [6]). For input 
domains where many sources or a high \z\ range are required, a different solver implementation may be 
needed, such as a varying precision arithmetic computation, working at higher floating-point accuracies. 

After sources' locations are recovered by taking the roots of the polynomials, the function 0^ of eq. 
^ is calculated, given the Newmann boundary conditions in eq. This was done using the Method 
of Fundamental Solutions (MFS), see e.g. [2^, where the solution to the Laplace for is approximated 
by a sum of (real) source terms, i.e. functions of the form cj){r) — ^ln|r — r^j, with the locations 
lying outside the domain D where 0l is harmonic. The charges Qi are then calculated to best fit the 
boundary conditions. 

Once the function is given everywhere in _D, the directions of the star-geodesics at each source 
are calculated, as described in Appendix A. The star-geodesics, emanating from the source were then 
traced by solving the differential equation, cq. ([2]), using a Runge-Kutta ODE solver. The cut-tree 
was constructed by taking one source at a time at some arbitrary order, and introducing a cut along a 
star-geodesic of the added source, that is closest to the cut-tree constructed so-far. This construction of 
the cut tree is admitantly arbitrary, and may not be optimal in some cases. 

Lastly, the algorithm described in section 0] is applied. 

6 Example Results 

In the first example, a domain with two boundary components is meshed, see fig. |6la. The side lengths 
of the outer and inner squares are 8 and 2.26 respectively, centered at the origin, and the inner square 
is rotated at angle a — O.SStt radians. The ratio of input resolution between outer and inner boundaries 
is 2.5. The inner and outer boundaries were sampled with 268 and 76 points, respectively, which were 
used both for calculating the Cauchy integral, and for the sources' locations calculation. Singularities of 
(j) are not required at the junctions points, since all inner angles are multiples of 7r/2 radians. As there 
is no flux of V(/) through either boundary (because the boundaries' curvature is k = 0, see eq. dH), and 
(j) is regular at the junctions), the total charge inside is zero (see remark [5]). Hence an equal number 
oi k — +1 and k — —1 charges were used. Fig. [SJa., shows the domain boundary, together with the 
calculated sources' locations. Squares mark sources with fc = -1-1, and circles sources with fc = — 1. The 
distribution is composed of 18 sources of each type, i.e. a total of 36 sources. The value of (p obtained 
by solving the Poisson equation with this charge distribution is shown in fig. [6lb., and compared to the 
input requested, denoted by ipinput- The difference A(/) = cj) — <j)input is also plotted. For visual clarity, 
the average of the ipmput was subtracted from both (p and (pinput (this is just a scaling of the resolution 
by a constant factor); this was also done in the resolution comparisons below, in fig. (I7]),c., fig. ([8]),b. 
and fig. (ini),c.. The number of sources was chosen to be the smallest for which A0 < 0.1 at all boundary 
points. This criterion is used for choosing the number of sources is also used in the following examples. 
The cut-tree used in creating the final mesh is shown in Fig. [Blc., and the final mesh in shown Fig. [Sid. 

The domain in the second example is the union of three unit radius circles which pass through the 
origin, see fig. [71a. The boundary has three junctions, all with the same inner angle, 9^^ — 57r/3. 
Boundary alignment requires (j) to be singular at the junctions, corresponding to sources with charge 
Q ~ — 7r/5, or fc-value k = Q/{t:/2) = —2/5, see eq. p7)) . (fT8]) . As explained in section [3.4.1 [ singular 
functions can be added the input resolution in order to avoid clustering of sources near the junction 
points. In this example, we compose the input </> from a sum of pairs of source terms, with fc-values as 
marked in fig. [7| ,a.. The pairs of nearby sources with opposite charge create the desired singularities at 
the junctions, and have only a small effect at distances large compared with the distance between the 
sources of each pair. 
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Figure 6: (Color online) Mesh generation procedure: a. Domain boundary (solid line) sampled at points 
(dots) as input to the IP algorithm. Resulting charge distribution is plotted: k = +1,-1 cone-points 
are plotted with squares and circles, respectively, b. </> obtained at boundaries. Points 1-268 correspond 
to the outer boundary component, the rest to the inner component. Input 0, obtained and their 
difference are plotted, c. Crosses at selected points (crosses), cut-tree composed af the boundary and 
star-geodesics (solid line), and selected additional star-geodesics (dashed lines), d. The final mesh. 
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We can now calculate the total flux of (f) through the boundary in this example. The (singular) flux 
through each junction is 7r/6 (positive value when V0 is directed outward). The flux through each of 
the three arcs is / dn<pds — — J dgOds — —Ad along the arc, and since for each arc A9 — 47r/3, the total 
flux through the boundary is 3 • (tt/G) + 3 • {—Att/3) — —7 (7r/2). This is in agreement with the formula 
given in remark [TJ Therefore, the sum of fc-values of the sources inside the domain should be —7, so 
n_ — n^- + 7. Hence, kin of eq. (|23p was —7; Zin = 1 was used. The lowest number of sources for which 
A(j) < 0.1, shown in fig. (Hb., has n+ — 2, = 9. The boundary was sampled at 525 points. 
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Figure 7: (Color online) a. Boundary with adjustment of resolution near junction points, b. Obtained 
source locations (blue), c. Input 0, obtained <j) and their difference, d. Final mesh. 



The boundary in the third example is composed of two boundary components, each with a single 
junction point. The inner angles at the inner and outer boundary junctions are 9in — 47r/3, 7r/2 respec- 
tively, so a singularity of <j) will form at the inner boundary junction. As in the previous example, we 
add a singular function to the input 0, by adding a pair of source terms, one at the junction point with 
k = 1/2, as obtained from eq. (fT7|) . (|18p with Om = 47r/3. No inner source term is required, and from 
the flux trough the hole boundary it follows that khoie — —3. After subtracting the junction source term 
with fc = 1/2 (see section a source with total charge of ky^oie — —3 — 1/2 is added inside the hole, 
see fig. Ela.. The inner and outer boundary components are sampled at 325 and 142 points, respectively. 
The distribution with = rt_ = 16, the smallest number of sources for which A<j) < 0.1, is plotted at 
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fig. [Hla., with squares and circles representing charges as in previous examples. Two pairs oi k = ±1 
charges are located inside the inner boundary, i.e., outside the domain. Their constitute a correction to 
the homogeneous (j)L, which is very small, since in each pair the opposite charges almost overlap. The 
obtained 4> reconstruction is given in fig. [Hlb.. The final mesh, fig. [SJc, contains 26 irregular vertices, 13 
of each type. This is because two pairs of sources lay outside the domain, and another pair (the pair of 
sources closest to the inner junction) was "eliminated" in the creation of the final mesh (section |4|) : the 
two opposite charges where united, giving a zero total charge. 
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Figure 8: (Color online) a. The boundary (solid line), sources to adjust input resolution near junc- 
tion (triangles), source distribution (squares and circles), and inside source location (red star), used in 
removing the harmonic part. b. Input ^, obtained </) and their difference, c. Final mesh. 



In the final example the boundary of the square [—2, 2] x [—2, 2] was assigned a varying input resolution 
proportional to: 

[y + 2)^ 

resolution oc 1.5 + (a; + 2)^ H . 

The resolution thus varies by a factor of 17 within the domain. As always, eq. (fTTj) . </> = In {resolution). 
The total charge is zero, so n^ — n_. The boundary was sampled at 536 points. A source distribution 
with Acf) < 0.1 was obtained with just 8 charges, i.e. with n^ = n_ — 4. Fig. [SJa. shows the source 
distribution. Fig. [51b. shows the obtained (p vs. the input (p, and the final mesh is shown in Fig. ^c. 
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It is instructive to compare the star-geodesics and cut-tree with the final mesh. The two are overlaid in 
Fig. inid. (this would create a "crowded" appearence when more sources are present). The black solid 
line traces the cut-tree, the dashed lines are additional star-geodesics, and the mesh is drawn with the 
background in gray lines. 
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Figure 9: (Color online) a. Boundary, source distribution (squares and circles), b. Input (f>, obtained (j) 
and their difference, c. Final mesh. d. Final mesh (gray lines), overlaid with the cut-tree (solid line) 
and additional star-geodesics (dashed lines). 



Cell shape quality was measured using a variant of the quality measure P |28j, as defined in |29| . 
(3 = 1 represents a square cell, while (3 = represents a cell with an inner angle of vr. Cells with a 
high aspect ratio are also given low /3-values. Table 1 shows the minimum and average (3 values, and 
the total number of cells, for the examples in Fig. For each example, the tabulated information is 
given for the mesh shown in the corresponding figure, and for a finer mesh of the same domain, where 
the resolution function was doubled (prodcing a mesh with about four times the number of cells) . Note 
that multiplying the resolution by a factor amounts to adding a constant to 0, which does not change 
the source distribution obtained from the IP algorithm, as this constant, which is part of h (z), is readily 
removed when h (z) is subtracted. 
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1062 


0.626 


0.964 




6059 


0.684 


0.987 



Table 1: Mesh statistics. 



7 Conclusions and Future work 

An unstructured quadrilateral mesh generation scheme in planar domains was presented. The method 
rests on a theoretical foundation, linking the mesh generation problem with the Inverse Poisson (IP) 
problem. An IP solution algorithm is presented, whose output is interpreted as the location and type 
(degree) of irregular vertices in the domain. The continuum fields obtained, describing mesh resolution 
and directionality, are conformal everywhere except on the irregular vertices, and fit the required input 
properties at the boundary, or at other user-defined locations. An algorithm for creating a valid final 
mesh is also presented. Example meshes feature irregular vertices where they are needed, in combination 
with highly regular regions where possible. 

Directions for future work include more sophisticated methods for solving the rational function inter- 
polation equations, and for constructing the final mesh. The relations between conformal unstructured 
mesh generation, the IP problem and rational function interpolation raise many research questions. These 
may lead to a deeper understanding of the properties of high quality meshes, and to better algorithms 
for creating them. 

Acknowledgements. The author would like to thank Mirela Ben-Chen, Shlomi Hillel, Dov Levine, 
Yair Shokef and Vincenzo Vitelli for helpful discussions and critical reading of the manuscript. 

8 Appendix A: Star-geodesic directions 

A cross is defined at every point p of the domain, that is not a singularity of the function 0. Geodesic 
curves that start at p can be drawn in all four directions of the cross. By the definition of a cross-field, 
such a geodesic will be aligned with the cross-field everywhere along the curve. If p is a singular point, 
with fc- value fcp 7^ 0, a cross is not defined at p, but there are 4 -I- fcp geodesies that are incident on p and 
follow the cross-field directions elsewhere on the curve. These will be called star- geodesies. For example, 
the geodesies drawn in fig. [6lc. and[9ld. are star-geodesics. 

Denote the angle from the x-axis around p by V', and the cross direction when p is approached from 
direction i/) by 6* (i/j), see fig. [101 To calculate the directions in which star-geodesics emanate from p, we 
first calculate the cross when p is approached from some direction, e.g. the positive x-axSs. This can be 
done using eq. ([3]) along a curve from the boundary to p which approaches p from the direction ■0 = 0. 
Once 9 (0) is known, 9 (ip) for any ^ can be calculated by using eq. Q along a small circular arc around 
p at radius r, a^. The singularity term at </) is | In |r|, and according to eq. ^ 

9{^)^9 (0) + / ^ds = (0) + ^^r = (0) + ^^r ^ 9 {Q) + (37) 
on or 4r 4 

Note that contributions to that are regular at p do not affect 9 (-0) for r ^ 0. 

The star-geodesic directions are those for which the cross is directed along the ray from p, or 

0(0)=0 + n|, (38) 



21 



Figure 10: Calculating the star-geodesics' directions. 

with n € Z. Substituting eq. ([37]) into eq ()38l ) and rearranging, we find 

'(O)-nf 



k 



1 4 



It is easy to show that there are exactly 4 + k different ip-values for which this equation is fulfilled. They 
are equally distributed around p. 
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